Exact results for quench dynamics and defect production in a two-dimensional model 



00 

o 
o 

(N 
C 

a 
in 



I 

C 

o 
o 



> 
(N 



o 



X 



K. Sengupta 1 , Diptiman Sen 2 and Shreyoshi Mondal 1 
1 TCMP division, Saha Institute of Nuclear Physics, 1/AF Bidhannagar, Kolkata 700 064, India 
2 Center for High Energy Physics, Indian Institute of Science, Bangalore, 560 012, India 

(Dated: February 2, 2008) 

We show that for a d-dimensional model in which a quench with a rate r _1 takes the system 
across a d — m dimensional critical surface, the defect density scales asti~ ]_/ T mv /(.* v + 1 ) > wnere v 
and z are the correlation length and dynamical critical exponents characterizing the critical surface. 
We explicitly demonstrate that the Kitaev model provides an example of such a scaling with d — 2 
and m = v = z = 1. We also provide the first example of an exact calculation of some multispin 
correlation functions for a two-dimensional model which can be used to determine the correlation 
between the defects. We suggest possible experiments to test our theory. 
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Quantum phase transitions have been studied exten- 
sively for several years Such a transition is accom- 
panied by diverging length and time scales [l[ leading to 
the absence of adiabaticity close to the quantum critical 
point. Thus the system fails to follow its instantaneous 
ground state when some parameter in its Hamiltonian is 
varied in time at a finite rate 1/r which takes the sys- 
tem across the critical point. Since the final state of the 
system does not conform to the ground state of its final 
Hamiltonian, defects are produced The defect den- 

sity n depends on the quench time rasn~ \ / T dv / i v z + x ) ^ 
where v and z are the correlation length and dynamical 
critical exponents at the critical point A theoretical 
study of such a quench dynamics requires a knowledge of 
the excited states of the system. Such studies have there- 
fore been mostly restricted to phase transitions in exactly 
solvable models in one or infinite dimensions 0, [(| 0, Hf • 
Experimental studies of defect production due to quench- 
ing of the magnetic field in a two-dimensional (2D) spin-1 
Bose condensate have been undertaken However, ex- 
act studies of quench dynamics have not been carried out 
so far for 2D spin models. 

In this Letter, we carry out such a study for the 2D 
Kitaev model [lfj . We show that when the quench takes 
the system across a critical (gapless) line, the density of 
defects scales as 1/ y/r. The Kitaev model has d = 2 and 
v = z = 1; hence, the above scaling is in contrast to 
the n ~ 1/t behavior expected when the system passes 
through a critical point [4[. In this context we provide 
a general discussion of the scaling of the defect density 
for a general <i-dimensional model with arbitrary v and 
z; we show that when the quench takes such a system 
through a d — m dimensional gapless (critical) surface, 
the defect density scales as n ~ \ j T mu / i zv + x ) , This re- 
sult is a generalization of the result of Ref. [3| and hence 
constitutes a significant extension of our current under- 
standing of defect production due to a quench. We also 
compute exactly some multispin correlation functions for 
our model, and use them to study the variation of defect 
correlations with the quench rate and the model param- 



eters. Such an exact analysis of defect correlations has 
not been carried out so far for 2D systems. 

The Kitaev model is a spin-1/2 model on a 2D honey- 
comb lattice with the Hamiltonian 



J2<T v j _ ll a v ]l + J 3 <tIi<tIi +1 ), 



j-\-l—even 



(i) 

where j and I denote the column and row indices of 
the lattice. This model, sketched in Fig. □] is known 
to have several interesting features [l2|, 13, fill It 
is a rare example of a 2D model which can be exactly 
solved 03, EL it supports a gapless phase for 
1^1 — J2 1 < 'h < J 1 + J 2 UfJ which is possibly con- 
nected to a spin liquid state and demonstrates fermion 
fractionalization at all energy scales [l3j]. In certain pa- 
rameter regimes, the ground state exhibits topological 
order and the low-energy excitations carry Abelian and 
non-Abelian fractional statistics; these excitations can 



be viewed as robust qubits in a quantum computer 11|. 
There have been proposals for experimentally realizing 
this model in systems of ultracold atoms and molecules 
trapped in optical lattices (l6| : such systems are known 
to provide easy access to the study of non-equilibrium 
dynamics of the underlying model. However, in spite of 
several studies of the phases and low-lying excitations of 
the Kitaev model, its non-equilibrium dynamics has not 
been studied so far. We will study what happens in this 
model when J3 is varied from —00 to 00 at a rate 1/r, 
keeping J\ and J2 fixed. 

One of the main properties of the Kitaev model which 
makes it theoretically attractive is that, even in 2D, it 
can be mapped onto a non-interacting fermionic model 
by a suitable Jordan- Wigner transformation 12|, |l4| , 



Hf = i 'Yl\ Jlbiia R-Mi + J2b ™ a n+M 2 + JaDfibnan], (2) 

ft 

where and are Majorana fermions sitting at the 
top and bottom sites respectively of a bond labeled n, 
n = ni + (2 i + §j) n 2 denote the midpoints of the 



2 




FIG. 1: Schematic representation of the Kitaev model on 
a honeycomb lattice showings the bonds Ji, J2 and J3. 
Schematic pictures of the ground states, which correspond 
to pairs of spins on vertical bonds locked parallel (antiparal- 
lel) to each other in the limit of large negative (positive) J3, 
are shown at one bond on the left (right) edge respectively. 
Mi and M2 are spanning vectors of the lattice, and a and b 
represent inequivalent sites. 
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FIG. 2: Plot of defect density n versus Jr and a = 
tan _1 (J2/Ji). The density of defects is maximum at Ji = J2. 



vertical bonds shown in Fig. [TJ and ni,ri2 run over all 
integers. The vectors ft form a triangular lattice. The 
x, y coordinates of the triangular lattice sites are given 
by x — ^/i{ni + n 2 /2) and y = 3n 2 /2. [We will refer 
to the sites of the honeycomb lattice either as (j, I) as 
in Eq. fT}, or as {a, ft) and (o, n) as in Eq. @; j + I 
is even (odd) for a (b) sites respectively.] The vectors 

Mi = + I j and M2 = 2 * — 1 3 are s P ann i n g vectors 
for the reciprocal lattice. The operator Dfi can take the 
values ±1 independently for each ft and commutes with 
Hp, so that the states can be labeled by the values of 
on each bond; t he g round state corresponds to = 1 
on all bonds [HE! EE EI irrespective of the sign of J3 
due to a special symmetry of the model [lj| • Since D,-% is 
a constant of motion, the dynamics of the model starting 
from the ground state never takes the system outside the 
manifold of states with D 7 i = 1 . 

For Dj{ = 1, Eq. (JSJ can be diagonalized as Hp = 
^k'^Hj^'ki where tpt = (at, bt) are Fourier trans- 
forms of an and bft, the sum over fc extends over half 
the Brillouin zone (BZ) of the triangular lattice formed 
by the vectors ft, and H? can be expressed in terms of 
the Pauli matrices a % (where er 3 is diagonal) as H% = 
2[ Ji sin(fc- Mi) - J 2 sm(k- M 2 )]ct 1 + 2[ J 3 + J x cos(fc • M x ) + 
J2 cos(/c ■ M2)]er 2 . The spectrum consists of two bands 

where 



with energies E~ 



Ej: = 2[{ Ji sin(fe • Mi 



J2sin(A:-M 2 )} 2 
+{ J 3 + Ji cos(k ■ Mi) + Ji cos(fc ■ M 2 )} 2 ] 1/2 (3) 



For \3\ — Ji\ < J3 < J1 + J2, the bands touch each other 
and the energy gap = Et 



values of k leading to a gapless phase 13, 12, 14, H)]. 

We will now quench Js(t) = Jt/r from —00 to 00 at a 
fixed rate 1/r, keeping J, J x and J 2 fixed at some positive 
values. The ground states of Hp corresponding to J3 — > 
—00(00), schematically shown in Fig. [TJ are gapped and 
have cj ; cr| l+1 = 1(— 1) for all lattice sites (j, I) of type b. 
To study the time evolution of the system, we note that 
after an unitary transformation U — exp(— ier 1 7r/4), we 
obtain Hp = ^j^^^HLipL, where HL = UH^U^ is given 
by Hi. = 2[Ji sin(fc ■ Mi) - J 2 sin(fc ■ Nh)]^ 1 + 2[J 3 (t) + 

J\ cos(fc • Mi) + J2 cos(fc • M 2 )]cr 3 . Hence the off-diagonal 
elements of Hi. remain time independent, and the quench 
dynamics reduces to a Landau-Zener problem for each 
k. The defect density can then be computed following a 
standard prescription 17 1: n = (1/-A) d 2 k p^, where 



E-. vanishes for special 



P£ = exp[-27rr{ Ji sin(A; • Mi) - J 2 sin(fc • M 2 )} 2 / J] (4) 

is the probability of defect production for the state la- 
beled by momentum fc, and A = 47r 2 /(3\/3) denotes the 
area of half the BZ over which the integration is carried 
out. A plot of n as a function of the quench time Jt 
and an angle a is shown in Fig. [2j here we have taken 
Ji[2] — J cos a [sin a]. We note that the density of defects 
produced is maximum when a = 7r/4 (Ji = J2). This oc- 
curs because the length of the gapless line through which 
the system passes during the quench is maximum for 
J\ = J2. Hence the system remains in the non-adiabatic 
state for the maximum time during the quench, leading 
to the maximum density of defects. Note that unlike 
some other models H, the defects here do not corre- 
spond to topological defects since the dynamics always 
keeps Df} = 1 on all bonds. 
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For sufficiently slow quench 2ir Jt 3> 1, pg is expo- 
nentially small for all values of k except near the line 
Ji sin(fc-Mi) = J2 sin(fc • M2); the contribution to the 
momentum integral in the expression for n comes from 
this region. Note that the line pg = 1 precisely corre- 
sponds to the zeros of the energy gap Ag as J3 is varied 
for fixed Ji, J2. By expanding pg about this line, we see 
that for a very slow quench, the defect density scales as 
n ~ 1/y/r. This demonstrates that the scaling of n with 
t crucially depends on the dimensionality of the critical 
surface since for a quench which takes the system through 
a critical point instead of a critical line, the defect de nsity 
of the Kitaev model, which has d = 2 and v = z = 1 [l8| , 
is expected to scale as 1/r [3| . This observation leads to 
the following general conclusion. 

Consider a d-dimensional model with v = z = 
1 which is described by a Hamiltonian Hd = 
Efc^ (cr 3 e(fc)t/r + A(fc)d+ + A*(k)o-) ip % , where a ± 

— (a 1 ±ia 2 )/2. Suppose that a quench takes the system 
through a critical surface of d — m dimensions. The de- 
fect density for a sufficiently slow quench is given by [TtJ 



=0.2 i . . 



. d d ke~ T ' Tf ^ 



- (1/40 I B z ddk ex Pi 

/ 2 , where Ad is the area of 
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n = (1/A d ) J BZ , 

the d-dimensional BZ, f(k) — |A(fc)| 2 /|e(fc)| vanishes 
on the d — m dimensional critical surface, a, [3 denote 
one of the m directions orthogonal to that surface, and 
g af 3 = [d 2 f(k)/dk a dkp\^ crmcalsmiacc . Note that this 
result depends only on the property that f(k) vanishes 
on a d — m dimensional surface, and not on the pre- 
cise form of f(k). For general values of v and z, we 
note that a Landau- Zener type of scaling argument yields 
A ~ l/ T zu /( zu + 1 ) ; where A is the energy gap When 
one crosses a d—m dimensional critical surface during the 
quench, the available phase space f2 for defect production 
scales as Cl ~ k m ~ A m/z - i/ T «W(^+i) ; this leads to 
n ~ \ I T mv I ^ zvJrV> . For a quench through a critical point 
where m — d, we retrieve the results of Ref. [4[ . 

Next we study the correlations between the defects 
produced during the quench. To this end, we define 
the operators Op = ibp^p^p. In the spin language, 
Og = For r ^ 0, Op can be written as a prod- 

uct of spin operators going from a b site at n to an a site 
at n + r: the product begins with a a x or a y at (j, I) and 
ends with a a x or a y at with a string of a z, s in 

between, where the forms of the initial and final a ma- 
trices depend on the positions of j + I and f + I'. Note 
that for J3 — ► —00(00), where the z component of each 
spin is locked with that of its vertically nearest neighbor, 

(^-00(00) |Of#-oo(oo)) = ±<W For tlie Kitaev model, 
it is known that the spin correlations between sites ly- 
ing on different bonds vanish, e.g., (o z a -cr^ - + -) = for 

13} . Therefore (Op) are the only non- vanishing 



^ 

two-point correlators of the model 1J| . A non-zero value 
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FIG. 3: Plot of (Op), sans the 5-function peak at the origin, 
as a function of r for four values of J2/J, for Jt — 5. 



the defects. In particular, a plot of (Op) versus r gives an 
estimate of the correlations between the defects. [Since 
Oi = 1, all the moments of Op can be found trivially: 
(Op 1 ) = (Op) if n is odd and = 1 if n is even.] 

After the quench, the system, for each momentum k, 
is described by a combination of tp_ OCl t with probability 
p^ and "000 £ with probability I ~ p^, where are 
the eigenstates of H? for J3 — > ±00. Hence (Op) can be 
computed in a straightforward manner: 



(O r ) 



5*3 



2 



<i 2 £; p£ cos(fc • r), (5) 



of (O^) for r ^ in the final state provides a signature of 



where the integral runs over half the BZ with area A. 

For large values of r, the dominant contribution comes 
from the region near the line J\ sin(fc • Mi) = J2 sin(fc • 
M2) where pt = 1. Introducing the variables fen and fcj^ 
which vary along and perpendicular to this line (along 
the directions ny and n± respectively), we see that the 
integrand in Eq. (O takes the form exp[— a(ko)rkj_ ± 
i(k + fc||7i|| + kj_hj_) ■ r\, where a(fc n ) is a number de- 
pending on fco- The evaluation of the integral over k± 
gives a factor of exp [— (r • njj 2 /(4a-r)] /y/ar. We thus 
see that the magnitude of the defect correlations goes 
as 1/y/r, while the spatial extent of the correlations 
goes as y/r. This is confirmed by the following relation: 
Er^(Or) = -2(V|pg)g =g = 247rr(J 2 + jf + J x J 2 )/J. 

To study the correlations between defects, we evaluate 
Eq. ([5]) numerically; the r dependence of (Op) is shown 
in Fig. [3] for several values of J2/J1, where J\ — J and 
Jt = 5. Here we have omitted the (^-function peak at 
r = in Eq. [5] so as to make the correlation at r ^ 
clearly visible. The plot of (Of) reflects the defect corre- 
lations. To understand the variation of these correlations 
with the ratio J2/J1 , we note that for large Jt, the max- 
imum contribution to (Op) comes from around the wave 



4 



0.08 



0.04 



0.00 



-0.04 



-0.08 




FIG. 4: Plot of (Of) at the points (1, 0) on the rn axis (black 
solid line), (0, 2) on the ri2 axis (blue dotted line), and (2, —2) 
along the —45° line in the ni — ri2 plane (red dashed line) as 
a function of a = tan _1 (J2/Ji), for J 2 = 1 and Jr = 5. 



vectors fco f° r which p(ko) = 1. For J 2 ^> (<C)Ji, this im- 
plies sm[k-M 2 (Mi)] = which yields k a ~ V3i± j. The 
maximum contribution to (Off) comes from cos(fco-r) = 1, 
i.e., k ■ r*= 0. Thus for J 2 » (<C)Ji, (Of) is expected 
to be maximal along the lines y = — (+)y/3x, namely, 
n\ = —n 2 (ni = 0) in the n\ — n 2 plane. This expectation 
is confirmed in Fig. [3] where (Off) can be seen to be max- 
imal along ri\ — —n 2 (ni = 0) for J 2 = 4(0.2)Ji. Such a 
strong anisotropy can be understood by noting that the 
Kitaev model reduces to a one-dimensional model when 
J 2 ^S> (<^.)Ji- For intermediate values of J1/J2, a grad- 
ual evolution of the defect correlations can be seen in 
Fig. [3] We note that if the Kitaev model can be real- 
ized using ultracold atoms in an optical lattice (ljj , such 
an evolution of the defect correlations with J1/J2 can, 
in principle, be experimentally detected by spatial noise 
correlation measurements as pointed out in Ref. [1 91 ] - 

We can obtain a different view of the spatial anisotropy 
of the defect correlations by studying (Off) as a function 
of a = tan _1 (J 2 />/i)- As a changes, the ratio J2/J1 
varies from to 00 while fixing J\ + j| = J 2 = 1 . A plot 
of (Off) at three representative points (m, n 2 ) = (1, 0) (on 
the ni axis), (0, 2) (on the n 2 axis), and (2, —2) (along the 
—45° line in the n\ — n 2 plane) as a function of a, shown 
in Fig. [U reveals the nature of the defect correlations. 
We see that as J 2 /Ji = tana is varied from to 00, the 
magnitude of the correlation at the point (1,0) on the 
n\ axis increases till it reaches a maximum at J\ = J 2 
(a = 7r/4), and then decays to as a approaches n/2. For 
the points (0, 2) on the n 2 axis and (2, —2) along the line 
with slope —45°, the correlation becomes maximum when 
J 2 <C Ji (a = 0) and J 2 >• J\ {a — tt/2) respectively, 



as expected from Fig. [3] We conclude that the spatial 
anisotropy of the correlations between the defects (Off) 
depends crucially on the ratio J 2 j J\. 

To conclude, we have shown that the density of defects 
produced during a quench of the Kitaev model through 
a critical line scales with the quench time as 1/y/r, in- 
stead of the 1/t behavior expected for a quench through 
a critical point. We have provided a general result for 
the defect density which reproduces the result of Ref. [H 
as a special case. We have also discussed the variation 
of the defect correlations with the model parameters and 
pointed out the possibility of detection of these varia- 
tions in experiments. These results significantly improve 
our general understanding of the scaling of the density of 
defects and their correlations in 2D systems. 

We thank A. Dutta and A. Polkovnikov for stimulating 
discussions. 
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